## Voting and racial attitudes (whites only) -- robustness checks for JOP
## Brian T. Hamel and Bryan Wilcox-Archuleta
## First: 11 March 2020  
## Last: 16 March 2020

## Loading packages
## install.packages(c("estimatr", "tidyverse", "magrittr", "texreg", "gridExtra", "scales"))
library(estimatr)
library(tidyverse)
library(magrittr)
library(texreg)
library(gridExtra)
library(scales)
library(lme4)

## Loading data
load("01_data/dta.RData")

## Number of people per zipcode
people_per_zip = dta %>%
        group_by(zipcode) %>%
        mutate(n = 1) %>%
        summarise(tot_people = sum(n, na.rm = TRUE)) 

mean(people_per_zip$tot_people)
sd(people_per_zip$tot_people)

people_per_zip %>%
        filter(tot_people >= 30)

## Re-estimating the main models with random slope and intercept
pres_dem = lmer(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density + (1 | zipcode),
                data = dta %>% filter(white == 1))

house_dem = lmer(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                + educ + pct_white + pct_black + pct_unemployed
                + pct_college + log_per_cap_inc + gini + south + non_rural
                + log_pop_density + (1 | zipcode),
                data = dta %>% filter(white == 1))

rr = lmer(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                + educ + pct_white + pct_black + pct_unemployed
                + pct_college + log_per_cap_inc + gini + south + non_rural
                + log_pop_density + (1 | zipcode),
                data = dta %>% filter(white == 1))

affirm = lmer(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                + educ + pct_white + pct_black + pct_unemployed
                + pct_college + log_per_cap_inc + gini + south + non_rural
                + log_pop_density + (1 | zipcode),
                data = dta %>% filter(white == 1))

## Table of coefs., and save
##############
## TABLE A4 ##
##############
texreg(list(pres_dem, house_dem, rr, affirm),
       file = "03_output/mlm.tex", 
       label = "mlm",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- Models with Random Intercept for Zip Code", 
       custom.model.names = c("\\textit{President}", "\\textit{U.S. House}", "\\textit{Racial Resentment}", 
                              "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1),
       include.loglik = FALSE,
       custom.gof.names = c(NA, NA, "\\# of Individuals", "\\# of Zip Codes", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.9)


## Adding past racial segregation -- 90
pres_dem_zseg90 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                              + educ + pct_white + pct_black + pct_unemployed
                              + pct_college + log_per_cap_inc + gini + south + non_rural
                              + log_pop_density + zipcode_dissim_90,
                              data = dta %>% filter(white == 1),
                              clusters = zipcode, se_type = "stata")

pres_dem_cseg90 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                              + educ + pct_white + pct_black + pct_unemployed
                              + pct_college + log_per_cap_inc + gini + south + non_rural
                              + log_pop_density + county_dissim_90,
                              data = dta %>% filter(white == 1),
                              clusters = zipcode, se_type = "stata")

house_dem_zseg90 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + zipcode_dissim_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

house_dem_cseg90 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + county_dissim_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

rr_zseg90 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + zipcode_dissim_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

rr_cseg90 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + county_dissim_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

affirm_zseg90 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + zipcode_dissim_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

affirm_cseg90 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + county_dissim_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

## Table of coefs., and save
##############
## TABLE A8 ##
##############
texreg(list(pres_dem_zseg90, pres_dem_cseg90,
            house_dem_zseg90, house_dem_cseg90,
            rr_zseg90, rr_cseg90,
            affirm_zseg90, affirm_cseg90),
       file = "03_output/seg90.tex", 
       label = "seg90",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- Controlling for Racial Segregation in 1990", 
       custom.model.names = c("\\textit{President}", "\\textit{President}",
                              "\\textit{U.S. House}", "\\textit{U.S. House}",
                              "\\textit{Racial Resentment}", "\\textit{Racial Resentment}",
                              "\\textit{Affirmative Action}", "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)",
                             "Zipcode Dissimilarity", "County Dissimilarity"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Adding past racial segregation -- 00
pres_dem_zseg00 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + zipcode_dissim_00,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

pres_dem_cseg00 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + county_dissim_00,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

house_dem_zseg00 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                             + educ + pct_white + pct_black + pct_unemployed
                             + pct_college + log_per_cap_inc + gini + south + non_rural
                             + log_pop_density + zipcode_dissim_00,
                             data = dta %>% filter(white == 1),
                             clusters = zipcode, se_type = "stata")

house_dem_cseg00 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                         + educ + pct_white + pct_black + pct_unemployed
                         + pct_college + log_per_cap_inc + gini + south + non_rural
                         + log_pop_density + county_dissim_00,
                         data = dta %>% filter(white == 1),
                         clusters = zipcode, se_type = "stata")

rr_zseg00 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density + zipcode_dissim_00,
                      data = dta %>% filter(white == 1),
                      clusters = zipcode, se_type = "stata")

rr_cseg00 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density + county_dissim_00,
                      data = dta %>% filter(white == 1),
                      clusters = zipcode, se_type = "stata")

affirm_zseg00 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                          + educ + pct_white + pct_black + pct_unemployed
                          + pct_college + log_per_cap_inc + gini + south + non_rural
                          + log_pop_density + zipcode_dissim_00,
                          data = dta %>% filter(white == 1),
                          clusters = zipcode, se_type = "stata")

affirm_cseg00 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                          + educ + pct_white + pct_black + pct_unemployed
                          + pct_college + log_per_cap_inc + gini + south + non_rural
                          + log_pop_density + county_dissim_00,
                          data = dta %>% filter(white == 1),
                          clusters = zipcode, se_type = "stata")

## Table of coefs., and save
##############
## TABLE A9 ##
##############
texreg(list(pres_dem_zseg00, pres_dem_cseg00,
            house_dem_zseg00, house_dem_cseg00,
            rr_zseg00, rr_cseg00,
            affirm_zseg00, affirm_cseg00),
       file = "03_output/seg00.tex", 
       label = "seg00",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- Controlling for Racial Segregation in 2000", 
       custom.model.names = c("\\textit{President}", "\\textit{President}",
                              "\\textit{U.S. House}", "\\textit{U.S. House}",
                              "\\textit{Racial Resentment}", "\\textit{Racial Resentment}",
                              "\\textit{Affirmative Action}", "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)",
                             "Zipcode Dissimilarity", "County Dissimilarity"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Adding past racial conflict
pres_dem_conflict = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density + goldwater + protest,
                     data = dta %>% filter(white == 1),
                     clusters = zipcode, se_type = "stata")

house_dem_conflict = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                              + educ + pct_white + pct_black + pct_unemployed
                              + pct_college + log_per_cap_inc + gini + south + non_rural
                              + log_pop_density + goldwater + protest,
                              data = dta %>% filter(white == 1),
                              clusters = zipcode, se_type = "stata")

rr_conflict = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                              + educ + pct_white + pct_black + pct_unemployed
                              + pct_college + log_per_cap_inc + gini + south + non_rural
                              + log_pop_density + goldwater + protest,
                              data = dta %>% filter(white == 1),
                              clusters = zipcode, se_type = "stata")

affirm_conflict = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                              + educ + pct_white + pct_black + pct_unemployed
                              + pct_college + log_per_cap_inc + gini + south + non_rural
                              + log_pop_density + goldwater + protest,
                              data = dta %>% filter(white == 1),
                              clusters = zipcode, se_type = "stata")

## Table of coefs., and save
##############
## TABLE A10 ##
##############
texreg(list(pres_dem_conflict, house_dem_conflict, rr_conflict, affirm_conflict),
       file = "03_output/conflict.tex", 
       label = "conflict",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- Controlling for Past Racial and Political Conflict", 
       custom.model.names = c("\\textit{President}", "\\textit{U.S. House}", "\\textit{Racial Resentment}", 
                              "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)",
                             "Support for Goldwater", "Civil Rights Protest"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Adding past racial income gap -- 90
pres_dem_zinc90 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + zipcode_inc_gap_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

pres_dem_cinc90 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + county_inc_gap_90,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

house_dem_zinc90 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                             + educ + pct_white + pct_black + pct_unemployed
                             + pct_college + log_per_cap_inc + gini + south + non_rural
                             + log_pop_density + zipcode_inc_gap_90,
                             data = dta %>% filter(white == 1),
                             clusters = zipcode, se_type = "stata")

house_dem_cinc90 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                         + educ + pct_white + pct_black + pct_unemployed
                         + pct_college + log_per_cap_inc + gini + south + non_rural
                         + log_pop_density + county_inc_gap_90,
                         data = dta %>% filter(white == 1),
                         clusters = zipcode, se_type = "stata")

rr_zinc90 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density + zipcode_inc_gap_90,
                      data = dta %>% filter(white == 1),
                      clusters = zipcode, se_type = "stata")

rr_cinc90 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density + county_inc_gap_90,
                      data = dta %>% filter(white == 1),
                      clusters = zipcode, se_type = "stata")

affirm_zinc90 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                          + educ + pct_white + pct_black + pct_unemployed
                          + pct_college + log_per_cap_inc + gini + south + non_rural
                          + log_pop_density + zipcode_inc_gap_90,
                          data = dta %>% filter(white == 1),
                          clusters = zipcode, se_type = "stata")

affirm_cinc90 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                          + educ + pct_white + pct_black + pct_unemployed
                          + pct_college + log_per_cap_inc + gini + south + non_rural
                          + log_pop_density + county_inc_gap_90,
                          data = dta %>% filter(white == 1),
                          clusters = zipcode, se_type = "stata")

## Table of coefs., and save
###############
## TABLE A11 ##
###############
texreg(list(pres_dem_zinc90, pres_dem_cinc90,
            house_dem_zinc90, house_dem_cinc90,
            rr_zinc90, rr_cinc90,
            affirm_zinc90, affirm_cinc90),
       file = "03_output/inc90.tex", 
       label = "inc90",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- Controlling for the Racial Income Gap in 1990", 
       custom.model.names = c("\\textit{President}", "\\textit{President}",
                              "\\textit{U.S. House}", "\\textit{U.S. House}",
                              "\\textit{Racial Resentment}", "\\textit{Racial Resentment}",
                              "\\textit{Affirmative Action}", "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)",
                             "Zipcode White-Black Income Gap", "County White-Black Income Gap"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Adding past racial income gap -- 00
pres_dem_zinc00 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + zipcode_inc_gap_00,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

pres_dem_cinc00 = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density + county_inc_gap_00,
                            data = dta %>% filter(white == 1),
                            clusters = zipcode, se_type = "stata")

house_dem_zinc00 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                             + educ + pct_white + pct_black + pct_unemployed
                             + pct_college + log_per_cap_inc + gini + south + non_rural
                             + log_pop_density + zipcode_inc_gap_00,
                             data = dta %>% filter(white == 1),
                             clusters = zipcode, se_type = "stata")

house_dem_cinc00 = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                         + educ + pct_white + pct_black + pct_unemployed
                         + pct_college + log_per_cap_inc + gini + south + non_rural
                         + log_pop_density + county_inc_gap_00,
                         data = dta %>% filter(white == 1),
                         clusters = zipcode, se_type = "stata")

rr_zinc00 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density + zipcode_inc_gap_00,
                      data = dta %>% filter(white == 1),
                      clusters = zipcode, se_type = "stata")

rr_cinc00 = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density + county_inc_gap_00,
                      data = dta %>% filter(white == 1),
                      clusters = zipcode, se_type = "stata")

affirm_zinc00 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                          + educ + pct_white + pct_black + pct_unemployed
                          + pct_college + log_per_cap_inc + gini + south + non_rural
                          + log_pop_density + zipcode_inc_gap_00,
                          data = dta %>% filter(white == 1),
                          clusters = zipcode, se_type = "stata")

affirm_cinc00 = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                          + educ + pct_white + pct_black + pct_unemployed
                          + pct_college + log_per_cap_inc + gini + south + non_rural
                          + log_pop_density + county_inc_gap_00,
                          data = dta %>% filter(white == 1),
                          clusters = zipcode, se_type = "stata")

## Table of coefs., and save
###############
## TABLE A12 ##
###############
texreg(list(pres_dem_zinc00, pres_dem_cinc00,
            house_dem_zinc00, house_dem_cinc00,
            rr_zinc00, rr_cinc00,
            affirm_zinc00, affirm_cinc00),
       file = "03_output/inc00.tex", 
       label = "inc00",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- Controlling for the Racial Income Gap in 2000", 
       custom.model.names = c("\\textit{President}", "\\textit{President}",
                              "\\textit{U.S. House}", "\\textit{U.S. House}",
                              "\\textit{Racial Resentment}", "\\textit{Racial Resentment}",
                              "\\textit{Affirmative Action}", "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)",
                             "Zipcode White-Black Income Gap", "County White-Black Income Gap"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Subsetting to above median % white
pres_dem_median = lm_robust(pres_dem ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density,
                     data = dta %>% filter(white == 1 & pct_white >= median(dta$pct_white, na.rm = TRUE)),
                     clusters = zipcode, se_type = "stata")

house_dem_median = lm_robust(house_dem ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density,
                            data = dta %>% filter(white == 1 & pct_white >= median(dta$pct_white, na.rm = TRUE)),
                            clusters = zipcode, se_type = "stata")

rr_median = lm_robust(mean_rr ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density,
                            data = dta %>% filter(white == 1 & pct_white >= median(dta$pct_white, na.rm = TRUE)),
                            clusters = zipcode, se_type = "stata")

affirm_median = lm_robust(affirm ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                            + educ + pct_white + pct_black + pct_unemployed
                            + pct_college + log_per_cap_inc + gini + south + non_rural
                            + log_pop_density,
                            data = dta %>% filter(white == 1 & pct_white >= median(dta$pct_white, na.rm = TRUE)),
                            clusters = zipcode, se_type = "stata")

## Table of coefs., and save
###############
## TABLE A13 ##
###############
texreg(list(pres_dem_median, house_dem_median, rr_median, affirm_median),
       file = "03_output/median.tex", 
       label = "median",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- % White > Median", 
       custom.model.names = c("\\textit{President}", "\\textit{U.S. House}", "\\textit{Racial Resentment}", 
                              "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "% Black Workers", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Subsetting to above 75th percentile % white
pres_dem_75 = lm_robust(pres_dem ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density,
                     data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.75, na.rm = TRUE)),
                     clusters = zipcode, se_type = "stata")

house_dem_75 = lm_robust(house_dem ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                        + educ + pct_white + pct_black + pct_unemployed
                        + pct_college + log_per_cap_inc + gini + south + non_rural
                        + log_pop_density,
                        data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.75, na.rm = TRUE)),
                        clusters = zipcode, se_type = "stata")

rr_75 = lm_robust(mean_rr ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                        + educ + pct_white + pct_black + pct_unemployed
                        + pct_college + log_per_cap_inc + gini + south + non_rural
                        + log_pop_density,
                        data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.75, na.rm = TRUE)),
                        clusters = zipcode, se_type = "stata")

affirm_75 = lm_robust(affirm ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                        + educ + pct_white + pct_black + pct_unemployed
                        + pct_college + log_per_cap_inc + gini + south + non_rural
                        + log_pop_density,
                        data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.75, na.rm = TRUE)),
                        clusters = zipcode, se_type = "stata")

## Table of coefs., and save
###############
## TABLE A14 ##
###############
texreg(list(pres_dem_75, house_dem_75, rr_75, affirm_75),
       file = "03_output/p75.tex", 
       label = "p75",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- % White > 75th Percentile", 
       custom.model.names = c("\\textit{President}", "\\textit{U.S. House}", "\\textit{Racial Resentment}", 
                              "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "% Black Workers", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Subsetting to above 90th percentile % white
pres_dem_90 = lm_robust(pres_dem ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density,
                     data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.9, na.rm = TRUE)),
                     clusters = zipcode, se_type = "stata")

house_dem_90 = lm_robust(pres_dem ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                        + educ + pct_white + pct_black + pct_unemployed
                        + pct_college + log_per_cap_inc + gini + south + non_rural
                        + log_pop_density,
                        data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.9, na.rm = TRUE)),
                        clusters = zipcode, se_type = "stata")

rr_90 = lm_robust(mean_rr ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                        + educ + pct_white + pct_black + pct_unemployed
                        + pct_college + log_per_cap_inc + gini + south + non_rural
                        + log_pop_density,
                        data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.9, na.rm = TRUE)),
                        clusters = zipcode, se_type = "stata")

affirm_90 = lm_robust(affirm ~ wac_pct_black + pid7 + ideo5 + female + age + faminc
                        + educ + pct_white + pct_black + pct_unemployed
                        + pct_college + log_per_cap_inc + gini + south + non_rural
                        + log_pop_density,
                        data = dta %>% filter(white == 1 & pct_white >= quantile(dta$pct_white, 0.9, na.rm = TRUE)),
                        clusters = zipcode, se_type = "stata")

## Table of coefs., and save
###############
## TABLE A15 ##
###############
texreg(list(pres_dem_90, house_dem_90, rr_90, affirm_90),
       file = "03_output/p90.tex", 
       label = "p90",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Whites) --- % White > 90th Percentile", 
       custom.model.names = c("\\textit{President}", "\\textit{U.S. House}", "\\textit{Racial Resentment}", 
                              "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "% Black Workers", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)